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This paper presents an integrated flight dynamic modeling method for flexible aircraft that captures cou- 
pled physics effects due to inertial forces, aeroelasticity, and propulsive forces that are normally present in 
flight. The present approach formulates the coupled flight dynamics using a structural dynamic modeling 
method that describes the elasticity of a flexible, twisted, swept wing using an equivalent beam-rod model. The 
structural dynamic model allows for three types of wing elastic motion: flapwise bending, chordwise bending, 
and torsion. Inertial force coupling with the wing elasticity is formulated to account for aircraft acceleration. 
The structural deflections create an effective aeroelastic angle of attack that affects the rigid-body motion of 
flexible aircraft. The aeroelastic effect contributes to aerodynamic damping forces that can influence aerody- 
namic stability. For wing-mounted engines, wing flexibility can cause the propulsive forces and moments to 
couple with the wing elastic motion. The integrated flight dynamics for a flexible aircraft are formulated by 
including generalized coordinate variables associated with the aeroelastic-propulsive forces and moments in 
the standard state-space form for six degree-of-freedom flight dynamics. A computational structural model 
for a generic transport aircraft has been created. The eigenvalue analysis is performed to compute aeroelas- 
tic frequencies and aerodynamic damping. The results will be used to construct an integrated flight dynamic 
model of a flexible generic transport aircraft. 


I. Introduction 

Modern aircraft are increasingly designed to be highly maneuverable in order to achieve high-performance mission 
objectives. Toward this goal, aircraft designers have been adopting light-weight, flexible, high aspect ratio wings in 
modern aircraft. Aircraft design concepts that take advantage of wing flexibility to increase maneuverability have been 
investigated. By twisting a wing structure, an aerodynamic moment can be generated to enable an aircraft to execute a 
maneuver in place of the use of traditional control surfaces. For example, a rolling moment can be induced by twisting 
the left and right wings in the opposite direction. Similarly, a pitching moment can be generated by twisting both 
wings in the same direction. Wing twisting or warping for flight control is not a new concept and was used in the 
Wright Flyer in the 1903. The U.S. Air Force conducted the Active Flexible Wing program in the 1980’s and 1990’s to 
explore potential use of leading edge slats and trailing edge flaps to increase control effectiveness of F-16 aircraft for 
high speed maneuvers. 1 In the recent years, the Active Aeroelastic Wing research program also investigated a similar 
technology to induce wing twist in order to improve roll maneuverability of F/A-18 aircraft. 2 

Structural deflections of lifting surfaces interact with aerodynamic forces to create aeroelastic coupling that can 
affect aircraft performance. Understanding these effects can improve the prediction of aircraft flight dynamics and can 
provide insight into how to design a flight control system that can reduce aeroelastic interactions with a rigid-body 
flight controller. Generally, high aspect ratio lifting surfaces undergo a greater degree of structural deflections than 
low aspect ratio lifting surfaces. In general, a wing section possesses a lower stiffness than a horizontal stabilizer or a 
vertical stabilizer. As a result, its natural frequency is normally present inside a flight control frequency bandwidth that 
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potentially can result in flight control interactions. For example, when a pilot commands a roll maneuver, the aileron 
deflections can cause one or more elastic modes of the wings to excite. The wing elastic modes can result in changes 
to the intended aerodynamics of the wings, thereby potentially causing undesired aircraft responses. Aeroservoelastic 
filtering is a traditional method for suppressing elastic modes, but this usually comes at an expense in terms of reducing 
the phase margin in a flight control system. 3 If the phase margin is significantly reduced, aircraft responses may 
become more sluggish to pilot commands. Consequently, with a phase lag in the control inputs, potential pilot- 
induced oscillations (PIOs) can occur. Numerous studies have been made to increase the understanding of the role of 
aeroservoelasticity in the design of flight control systems. 3-6, 15 

Aeroelasticity remains a fertile field of research in spite of many early contributions to this field dated back to the 
1930’s. In the early days of National Advisory Committee for Aeronautics (NACA), numerous studies on aeroelas- 
ticity by investigators, such as Theodorsen and Garrick, had made significant advances in this field. 7-11 Theodorsen’s 
theory of unsteady aerodynamics still remains an essential tool for aeroelastic analysis. Recently, interests in coupled 
flight-structural dynamics have seen a renewal. Due to flexibility of modern aircraft structures, flight dynamic models 
of rigid-body aircraft have limitations and cannot accurately predict behaviors of flexible aircraft when elastic modes 
participate in the rigid-body motion. Recently, some investigators have investigated theoretical approaches to devel- 
oping integrated flight dynamics with aeroelasticity. Shearer develops an integrated flight dynamic model for a rep- 
resentative High Altitude Long Endurance (HALE) vehicle. 13 Meirovitch and Tuzcu also develop another integrated 
approach to flight dynamics of flexible aircraft. 14, 15 In both of these approaches, detail inertial-aeroelastic-propulsive 
force coupling in the governing structural dynamic partial differential equations were not made available. 

The purpose of this study is to produce a modeling capability for integrated flight dynamics of flexible aircraft 
that can better predict some of the complex behaviors in flight due to multi-physics coupling. Some of the important 
features in the present method are: inertial force coupling due to aircraft rigid-body acceleration, bending-torsion 
coupling due to wing pre-wist, aeroelastic -propulsive force coupling due to engine mounting on flexible wing struc- 
tures, and lastly an extension of Theodorsen’s theory to include chordwise bending in the aeroelastic angle of attack. 
Through the use of generalized coordinates that represent elastic deflections, the standard flight dynamic equations for 
six degree-of-freedom motion are modified to include effects of aeroelasticity and propulsive forces on flexible wings. 

II. Reference Frames 




Fig. 1 - Aircraft Reference Frames 

Figure 1 illustrates three orthogonal views of a typical aircraft. Several reference frames are introduced to facilitate 
the rigid-body dynamic and structural dynamic analysis of the lifting surfaces. For example, the aircraft inertial 
reference frame A is defined by unit vectors ai, a2, and a3 fixed to the non-rotating earth. The aircraft body-fixed 
reference frame B is defined by unit vectors bi, b2, and b3. The reference frames A and B are related by three 
successive rotations: 1) the first rotation about a 3 by the heading angle i f/ that results in an intermediate reference 
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frame A defined by unit vectors a ; , a,, and a 3 (not shown), 2) the second rotation about a 0 by the pitch angle 9 
that results in an intermediate reference frame B defined by unit vectors , b 2 , and b, (not shown), and 3) the third 
rotation about bj by the bank angle 0 that results in the reference frame B. This relationship can be expressed as 
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The left wing elastic reference frame D is defined by unit vectors di, do, and d3. The reference frames B and 
D are related by three successive rotations: 1) the first rotation about br by the elastic axis sweep angle — A that 
results in an intermediate reference frame B defined by unit vectors b j , b^ and b 3 (not shown), 2) the second rotation 
about negative b , by the elastic axis dihedral angle T that results in an intermediate reference frame D defined by unit 
vectors d,, d 2 , and d 3 (not shown), and 3) the third rotation about dj by an angle n that results in the reference frame 
D. This relationship can be expressed as 
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Generally, the effect of the dihedral angle can be significant. A full analysis with the dihedral angle can be 
performed but can also result in a very complex analytical formulation. Thus, to simplify the analysis, the dihedral 
effect is assumed to be negligible in this study. The right wing reference frame C can be established in a similar 
manner. In the analysis, the aeroelastic effects on the fuselage, horizontal stabilizers, and vertical stabilizer are not 
considered, but the analytical method can be formulated for analyzing these lifting surfaces if necessary. In general, a 
whole aircraft analysis approach should be conducted to provide a comprehensive assessment of the effect of flexibility 
on aircraft stability. Such an analysis is feasible in the current framework as described in this study. 


III. Elastic Analysis 

In the subsequent analysis, the combined motion of the left wing is considered. The wing has a varying pre- 
twist angle y(x) common in many aircraft. Typically, the wing pre-twist angle varies from being nose-up at the wing 
root to nose-down at the wing tip. The nose-down pre-twist at the wing tip is designed to delay stall onsets. Under 
aerodynamic forces and moments, wing structural deflections introduce strains in the wing structure. For high aspect 
ratio wings, an equivalent beam approach can be used to analyze structural deflections with a reasonable accuracy. 
The equivalent beam approach is a typical formulation in many aeroelasticity studies. 12,16 Experimental validation 
can show that equivalent beam approach is accurate for an aspect ratio as low as 3:1. The internal structure of a wing 
typically comprises a complex arrangement of load carrying spars and wing boxes. Nonetheless, the elastic behavior 
of a wing can be captured by the use of equivalent stiffness properties. These properties can be derived from structural 
certification testing that yields information about wing deflection as a function of loading. 

Consider an airfoil section on the left wing as shown in Fig. 2 undergoing bending and twist deflections. 
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Fig. 2 - Left Wing Reference Frame 
Let (x,y,z) be the coordinates of a point Q on the airfoil. Then 
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where t] and ^ are local airfoil coordinates, and 7 is the wing section pre-twist angle, positive nose-down. 

The axial or extensional deflection of a wing is generally very small and therefore can usually be neglected. Let V 
and W be chordwise and flapwise bending displacements of point Q, respectively, and let © be a torsional twist angle 
about the x-axis, positive nose-down. Then, the rotation angle due to the structural deformation can be expressed as 


8 (x, t ) = ©d 1 - IV, -d 2 + V x d 3 


( 4 ) 


where the subscripts x and t denote the partial derivatives of V , IV, and ©. 

Let (x 1 , y 1 , z 1 j be the coordinates of point Q on the airfoil in the reference frame D. Then the coordinates (x 1 . y 1 , z 1 ) 
are computed using the small angle approximation as 11 
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Differentiating xi, y 1, and zi with respect to x yields 
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Neglecting the transverse shear effect, the longitudinal strain is computed as 
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Therefore, the longitudinal strain is expressed as 


_ si£_ _ -yVxx ~ zW xx + ( y 2 + z 2 ) y®x 
s * 1 + \ (y 2 + z 2 ) (7) 2 

For a small wing twist angle y, the longitudinal strain is further simplified as 
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l» 2 + z ! ) (y+0,) 


(ID 


( 12 ) 


’ M x ' 


GJ&x 



My 

= 

0 


[ Ee 

. M z . 


0 


J 


— z 
-y 


dydz = E 


-B 2 y 

-B 37 


-B 2 r 

-ByJ 

Iyy 

~hz 

~Iyz 

hz 



' ©, " 


w„ 


. . 


where the area moments of inertia and the bending-torsion coupling constants B \ , B 2 , and By are defined as 
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The strain analysis shows that for a pre-twisted wing the flexural deflections V and W are coupled to the torsional 
deflection © via the slope of the wing pre-twist angle. This coupling is usually not present in the Euler-Bernoulli beam 
formulation, and can be significant if the term 7 is dominant as in highly twisted wings such as turbomachinery rotor 
blades. 


IV. Inertial Coupling 

During a high-g maneuver, rigid-body acceleration of an aircraft due to the Coriolis effect can generate inertial 
forces on lifting surfaces. For highly flexible wings, the inertial forces can couple with the structural deflection. 
In particular, if an aircraft experiences a very rapid roll rate, the inertial forces created by the roll acceleration can 
cause natural frequencies to change due to the stiffness contribution of the roll acceleration. This is a well-known 
phenomenon for rotary wing structures, whereby the rotation contributes to the apparent structural stiffness which is 
known as rotational stiffening. Generally, the roll motion of an aircraft has a greater rotational stiffening effect than 
the pitch and yaw motions due to a larger moment arm, which is the wing span. 

For the analysis, only the roll motion of an aircraft is considered. The kinematic relationship is to be developed 
to establish the velocity and acceleration due to rigid-body aircraft motion at a reference point on the wing. The 
velocity establishes aeroelastic angle of attack that is used to develop aeroelastic forces and moments. Similarly, the 
acceleration establishes the inertial force coupling with the wing elasticity. 

A. Kinematics 

Working in the left wing reference frame D, we define points O, P, and Q as the center of gravity of the aircraft, the 
origin of the reference frame D of the left wing, and the reference point of a mass particle on the airfoil section as 
shown in Fig. 2. Then the position vectors from O to P and from P to Q are defined as 

r ° P = r° p b l +r° p b 2 + r° r b 2 =x P A x +y/>d 2 +Z/>d 3 (16) 

v PQ =x x (x,f)d! + 7 ! (x,f)d 2 +zi (x,t) d 3 (17) 
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xp —sin A —cos A 0 r* 

y P = —cos A sinA 0 r y p (18) 

_ zp \ [ 0 0 - 1 J [ r z P . 

Using Kane’s dynamical notations, the velocity of the center of gravity, point O, and the angular velocity of the 
aircraft are defined as 
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where u and w are the axial and vertical velocity components at point O , respectively; p is the roll rate; and (o x and G)y 
are the angular rates in the reference frame D 
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The acceleration of the center of gravity and angular acceleration of the aircraft are computed as 
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where ay, a y , and a? are the acceleration components at point O in the reference frame D 
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The velocity of point Q is contributed by the structural deflection and is computed as 
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The acceleration of point Q is also contributed by the structural deflection and is computed as 
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and a x , a y , and a z are the structural-deflection induced acceleration components in the reference frame D 
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It is observed that angular rates and accelerations are coupled with the structural deflection to give rise to additional 
acceleration components on a wing structure. 
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B. Inertial Forces and Moments 

For simplicity, we assume that the neutral axis coincides with the elastic axis. This is a reasonable approximation 
and good results have been obtained with this assumption. When the neutral axis is offset from the elastic axis, the 
expressions for inertial forces and moments can be quite complex. The offset between the neutral axis and elastic axis 
is evaluated as 
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The inertial forces at the elastic center due to the structural deflection per unit length are then computed as 
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where a*, a* and a* are the rigid-body acceleration components at the elastic center 
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The inertial moments due to the structural deflection per unit length are now computed as 
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V. Aeroelastic Analysis 

The relative velocity of the air approaching a wing section includes the contribution from the wing structural 
deflection that results in changes in the local angle of attack. Since aerodynamic forces and moments are dependent 
on the local angle of attack, the wing structural deflection will generate additional elastic forces and moments. The 
local angle of attack depends on the relative approaching air velocity as well as the rotation angle (? from Eq. (4). 
The relative air velocity in turn also depends on a structural-deflection induced velocity. The oscillation of the airfoil 
results in an unsteady circulation. Kussner showed that the velocity at the 3 /4-chord point determines the circulation 
force on an oscillating airfoil. 17 This is in contrast to the lift due to circulation that acts at the 1 /4-chord point in a 
steady motion. 


-b/2 , b/2 



AC = Aerodynamic Center 
EA = Elastic Axis 


Fig. 3 - Airfoil Coordinates 

Referring to Fig. 3, we adopt the convention used in Theodorsen’s theory of unsteady aerodynamics whereby the 
chord of an airfoil has a length 2b and the elastic center is located at a distance —ab from the mid-chord, where the 
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parameter a is between -1 and 1 and a is negative when the elastic center is forward of the mid-chord, which is usually 
the case for a conventional airfoil 1617 . The local velocity is computed by evaluating the velocity of point Q in Eq. 
(25) with the coordinates at the 3/4-chord point yc = (f — a) Z?cos y and zc = (j — a)/? sin yin the reference frame D 


(34) 


In order to compute the aeroelastic forces and moments, the velocity must be transformed from the reference 
frame D to the airfoil local coordinate reference frame defined by (ju , rj , | ) (see Fig. 2). Then the transformation can 
be performed using three successive rotation matrix multiplication operations as 
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where vjj, vjj, and v| are rigid-body velocity components given by 
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and Avjj, Av/ , and Av| are the induced velocity components due to the structural deflection given by 
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A. Aeroelastic Angle of Attack 


Referring to Fig. 4, the local aeroelastic angle of attack on the airfoil section is due to the velocity components and 
i'/ and is computed as 
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Fig. 4 - Aeroelastic Angle of Attack 


Fet a* be a local rigid-body angle of attack. Then a* is computed as 
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where we recognize that a is the aircraft angle of attack which is defined as 
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B. Aeroelastic Forces and Pitching Moment 

In unsteady aerodynamics, the lift force is comprised of noncirculatory and circulatory components. The non circu- 
latory component is due to the apparent mass and inertia effects which are generated when the wing has a non-zero 
acceleration. This acceleration causes the surrounding air which has a finite mass to generate inertial forces that oppose 
the acceleration. The circulatory component is more important for wing sections. This is due to the vortical strength 
of circulation that generates lift. In unsteady aerodynamics, the vortices are shed in the flow in a complex fashion. 
Unsteady thin-airfoil theories developed by Theodorsen and Peters can be used to estimate the effect of unsteady shed 
vorticity. The former theory is appropriate for classical flutter analysis while the latter is a finite-state theory cast in 
the time-domain and is appropriate for time-domain analysis. 18 

The total aeroelastic forces and pitch momenting at the elastic center are the sum of the noncirculatory and circu- 
latory components and are expressed as 

sn sc 

J y J y 

n + ft 

n c 

m x m x 

where the superscripts a, n, and c denote aeroelastic, noncirculatory, and circulatory, respectively. 




1. Noncirculatory Forces and Pitching Moment 

The noncirculatory unsteady forces are due to the apparent air mass effect which is based on the acceleration of the air 
mass enclosed by a circular cylinder whose diameter is the airfoil chord. These noncirculatory forces are computed 
from the acceleration in the reference frame D acting at the mid-chord as 
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where y = —a/? cosy, z = —absiny, and p TO is the air density. 

In addition, the structural rotation of the wing induces a downwash component 


Av| ~ —ucosAa e (47) 

where u is the free-stream air speed, which generates a noncirculatory force acting at the 3 /4-chord point equal to 


9 ^ Ar) d (Xt 

/g = -pooTtb 2 ^ ~ poo7tb~ucosA 


dt 


dt 


(48) 


The angular acceleration of the wing twist also induces a nose-down pitching moment at the mid-chord 



d 2 a e 
dt 2 


Hence, the total noncirculatory forces and pitching moment at the elastic center are computed as 


fn 

Jy 


/;-/|siny 

fz 

= 

cosy 

. < . 


_ yf"-zf" + (k~a)bfl+m n x _ 


If p = 0 and y is small, the noncirculatory forces and pitching moment are approximately equal to 


fy 


V„ T yab@tt 4- yu cos A ~yf~ 

fz 

= -pooTtb 2 

W tt — ab& tt — u cosA^jp 

_ m” _ 


—ab (W tt — ab&tt — yVtt ) — (| — a) bucosA—jf- — \b 2 f^- 


Equation (51) is in agreement with the well-established results in aeroelasticity. 17 


(49) 


(50) 


(51) 


2. Circulatory Aeroelastic Forces and Pitching Moment 

The circulatory lift, drag, and pitching moment based on the Theodorsen’s theory are given by 


d 

dx 
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\ 

C L 


i 

Of 

1 

D 

= Pocbu 2 cos 2 A < 

1 

Cd 

H - C (^) CC e 

Co, a 
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i 
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(52) 


’ Cl " 


Cl,0 + CL,a a * + C[g a 

Co 

— 

Cd, o + Co,aOC* + Co,s a &a + KCy 

Cm 


Cm,AC T C m fi n 8a 


The 2-D section lift, drag, and pitching moment coefficients Cl, Co, and C m due to the rigid-body angle of attack 
are defined as 


(53) 


where 8 a is the aileron deflection, K is the induced drag constant, and C, n / \c is the quarter-chord pitching moment 
coefficient which is positive nose up and independent of the angle of attack. 

The function C (k) is called the Theodorsen’s function which is a complex-valued function of the reduced frequency 
parameter k 

bco 

k= (54) 

u cos A 

where co is the flutter or aeroelastic mode frequency. 

( 2 ) 

C (k) can also be expressed in terms of Hankel functions of the second kind H„ (k) as 


C(k) = 


Hf\k) 


H 


( 2 ), 


■i '(k)+iH®(k) 

where F (k) > 0 and G (k) > 0. 

When k = 0, the airfoil motion is steady and C (k) is real and unity. As k increases, there is a phase lag introduced 
as the magnitude of G ( k ) increases as shown in Figs. 4 and 5. The limiting values of F ( k ) and G (k) are 1 / 2 and 0 as 
k^ 00.12 


= F(k) — iG{k) 


(55) 
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Fig. 5 - Theodorsen’s Function 



Fig. 6 - Theodorsen’s Function Phase lag 


The aerodynamic forces and moment on a 2-D wing section acting at the elastic center in the reference frame D 
are determined by 


' fy ' 


—L sin ( a c + 7 ) + D cos (a c + 7 ) 

n 

LI 

dx 

Lcos (a c + 7 ) +£>sin(a c + 7 ) 

. m x . 

— M — + a) b [Lcos a c + D sin a c ] 


(56) 


as illustrated in Fig. 7. 



Fig. 7 - Airfoil Forces and Moment 
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For a small elastic angle of attack, the circulatory unsteady aerodynamic forces and moment are approximated as 


fy 

n 

.,,C 


= poobu 2 cos 2 A 
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1 

8 

1 

1 

Q 

H- C (&) cc e 

C z ,cc 

\ 

2 bC x 


2 bC x .a 

J 


where C y , C z , and C x are the force and moment coefficients due to rigid-body aerodynamics 


Cy 


-Cl sin (a* + y) + C D cos (a* + y) 

Q 

= 

Cl cos (a* + y) + C D sin (a* + y) 

C x 


-C m ~ {3 + f) [Cl cos a* 

1 

* 

.5 

’c /2 

+ 


(57) 


(58) 


and Cy^a, C z , a , and C x ,a are the derivatives of the force and moment coefficients with respect to the angle of attack 


1 

s 

1 


— {Cl, a + C D ) sin (a* + y) + {Cp,a — Cl) cos ( 0 :* + y) 

Cz,a 

= 

{Cl, a + Co) cos (cC + y) + {Co.a — Cl ) sin (cx* + y) 


C x ,a 


1 

1 

+ 

( c L,a + c *d) cos a * + {Cd,cc ~ Cl) sin a* 



C. Aerodynamic Damping 

Dynamic stability of a system is dependent on dissipative forces acting on it. The dissipative forces contribute posi- 
tively to damping of the system. Dynamically stable systems require positive damping. Aeroelastic forces give rise to 
the aerodynamic damping mechanism that influences aerodynamic stability of an oscillating wing or a flight vehicle. 
Positive aerodynamic damping results in aerodynamically stable operation. Conversely, negative aerodynamic damp- 
ing causes flutters and self-excited vibrations. Flutter boundaries are defined by air speed at which the aerodynamic 
damping crosses from a positive value to a negative value. 

Consider the case when p = 0 and y = 0, the damping forces are obtained from the noncirculatory and circulatory 
components of the aeroelastic forces as 


fd 

Jy 

T 9 ( aVwtanA _ \ 

- 7 ry 

9 9 , \ 

1 

■P 

a 

1 

ft 

= Poo^ 2 wcosA( — - — VV^tanA — 0 r I 

\ cosA J 

K 

+ poobu" cos" AC (k) a e 

C z ,a 

. m i . 
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1 
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2 bC x ,a 


a J 

(60) 


For a harmonic motion, the structural deflections may be expressed as 


’ © ’ 


4>© (x) 

W 

= 

(*) 

V 




(61) 


where <t>©, <t>w, and <t>y are the elastic mode shapes of the torsion, flapwise bending, and chordwise bending modes, 
respectively. 

Then, the elastic angle of attack may be expressed as 


a e = ( a r + iaj)e ,cot 

where a r and a, are the real and imaginary parts of a e 


a r 


a<i> v tan A 
cosA 


— <J> w tanA — <f>© 


a,- = -co 


+ (7 - a) M>© - a*<Fy 
mcosA 


(62) 


(63) 

(64) 


where the prime denotes derivative with respect to x. 
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The damping forces and pitching moment are then obtained as the imaginary part of the circulatory aeroelastic 
forces and pitching moment 


ft 


— Jty 
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u* 
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ft 

= poob 2 ucosAa r a>ie ,0,t 
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+ poobu 2 cos 2 A [— G (k) a r + F (k) a,-] ie im 
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which can also be written as 


(65) 
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u* 

1 

+ Poobu cos AF (k) 
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( 66 ) 


It should be noted that the sign of the aerodynamic damping is influenced by the sign of Cl.u + Cd. While Cd 
is always positive, Cl,u can be negative when the airfoil is stalled. Therefore, it is possible that the aerodynamic 
damping can become negative. If the overall damping which includes structural damping that inherently exists in 
the wing structure transitions from a positive value to a negative value, the ensuing motion will be aerodynamically 
unstable due to positive work inputs to the wing by the air. For example, the generalized damping coefficient for the 
flapwise bending, by neglecting the term Cd.u — Cl and the contributions from chordwise bending and torsion, may 
be estimated as 

PocMcosAF (k) f b (Ci )Ce + Cd) cos a* (x) [4>w (x)] 2 t/x 


C w — 


(67) 


2 co f pA [<t>w (x)} 2 dx 

It is obvious that if C/ a becomes negative then if it possible for to be negative. The aerodynamic damping for 
the torsion is more complex. In general, the aerodynamic damping for the combined bending-torsion motion must be 
analyzed by matrix analysis 


VI. Gravity and Propulsive Forces and Moments 

The gravity and propulsive forces are significant contributing active forces acting on the aircraft airframe. These 
forces can influence structural deflections of a wing. For high-aspect ratio, flexible wings, gravity can significantly 
offset wing bending deflections at low airspeed. Similarly, the engine thrust force can also affect twist and bending 
deflections of flexible wings. 


A. Propulsive Forces and Moments 


To deal with the propulsive force, we assume that the left engine produces a thrust vector aligned with the bi -direction 
at the thrust center E in the rigid-body aircraft reference frame B. Point E is assumed to be located relative to point P 
at x = Xe, y = yE, and z = Ze forward and below the elastic center of the wing section in the reference frame D. Since 
the engine is mounted on a flexible wing, then the thrust center E is dependent on the wing structural deflections. 

Since the engine thrust is a concentrated force, it can be formulated as a distributed force using the Dirac delta 
function which is defined as 


5(x-x 0 ) 


00 ,x = xo 
0 ,x^x 0 


( 68 ) 


such that for an arbitrary function g (x), then 


J g(x) 8 (x-xo)dx = g (xo) 


( 69 ) 
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Then, the distributed forces in the reference frame D due to the small structural rotation angle 8 are computed 
using three successive rotation matrix multiplication operations as 
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— cos A 
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— cos A — V v sinA 
— W.sin A — ©cosA 


(70) 


where T is the engine thrust. 

The distributed thrust forces are transformed back into the aircraft reference frame B as 
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— sin A 

— cos A 
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— sinA + V v cosA 
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« 8 (x — xe) T 

— cos A 

sin A 

0 


— cos A — V x sin A 

« 8 (x — xe) T 

-v x 
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0 

0 

-1 


— ITvsinA — ©cosA 


W T sinA + ©cosA 


The left engine thrust force can now be obtained by integration as 

T L = J 8 (x — xe) T [bi — V x b 2 — (VV^sin A + ©cos Aj ba] dx = T [bi — V x b 2 + (W v sinA + ©cosA)b3] A . =V£ 


(71) 


(72) 


It can be observed that the structural deflections at the thrust center generate additional thrust force components 
in the bo- and b 3 -directions. Thus, the coupled effect between the propulsive force and structural deflections can be 
significant if the wing is highly flexible and the engine thrust is large. 

The distributed moments due to the propulsive force in the reference frame D are computed as 


(73) 


B. Gravity Forces and Moments 

The gravity can exert significant forces and moments on a wing. The gravity forces include the weight of the wing 
including fuel, and the engine weight. The distributed gravity forces can be expressed in the reference frame D as 
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(74) 


where p is the mass density of the wing including fuel, A is the cross sectional area, and ihe is the engine mass. 

Assuming that the center of gravity of the wing and fuel coincides with the elastic axis, and the center of gravity 
of the engine coincides with the thrust center, then the distributed moments due to gravity in the reference frame D are 
computed as 


(75) 


m s x 


—yE cos 0 cos 0 — ze (cos A + cos 0 sin 0 sin A) 

mf, 

= 8 (x-XE)niEg 

— ze (sin 0 sin A — cos 0 sin 0 cos A) 

. m t . 


—yE (sin 0 sin A — cos 0 sin 0 cosA) 


VII. Coupled Structural Dynamic Equations 


The equilibrium conditions describe the force and moment balance of all forces and moments acting on a wing. 
The resulting force and moment equilibrium conditions are given by 


' F x ' 
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d_ 

dx 


M x ' 


m x 


’ F z V x -F y W x ' 


' 0 " 

My 

+ 

m y 

+ 

—F X W X + F Z 

= 

0 
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. m z . 


_ - F x V x + F y _ 


_ 0 _ 


( 77 ) 


where f x , f y , f z , m x , m y , and m z are distributed forces and moments due to inertial, aeroelastic, gravity, and propulsive 
effects 

f(x,y,z) = f(x,y,z) + f(x,y,z) + f(x,; y,z) + ^(w) (78 ) 

m ( W ) = m \x,y,z) + m (x,y,z) + m (x,y,z) + m (w) ( 79 ) 

The shear forces F y and F z can be eliminated by solving the last two equations in Eq. (77). Neglecting nonlinear 
terms, the equilibrium conditions can be written as 


dM x 

dx 


+ m x - m*V x + m*W x = 0 


(80) 


d 2 M y 

dx 2 

d 2 M z 

dx 2 



(81) 

(82) 


where the superscript * denotes rigid-body forces and moments, and a is a dummy variable that replaces x. 

The resulting equations are three structural dynamic partial differential equations that relate the flapwise bending, 
chordwise bending, and torsion with the coupled effects under consideration. The integral term involving the rigid- 
body inertial force / A *can be important if the inertial acceleration of the aircraft is significant. This term gives rise to 
the rotational stiffening effect when the angular speed of a wing structure is large. 


A. Example 

Consider a flight vehicle with non-twisted, unswept wings y = 0, A = 0, on a horizontal flight, 0=0, making a 
constant roll rate p maneuver. The rigid-body and elastic angles of attack from Eqs. (39) and (42) reduce to 


a = a 


- p(x P +x ) 


(83) 


a _ _w r £^_ 0 _. Wf + (2- fl ) fc0 ' 


(84) 


u u cos A 

Neglecting the chordwise bending motion and the damping forces, then the structural dynamic equations for flap- 
wise bending and torsion are 
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For a quasi-steady-state motion when the wing responds statically to the inertial, aeroelastic, and propulsive forces, 
by setting all the partial time derivative terms to zero in the above equations, one obtains 


— (GJ& X ) X 


PP Izz® 


pocb 2 u 2 cos 2 AC Zi a 



poobu 2 cos 2 KlbC x 


(87) 


{EIyyW xx ) xx - pAp 2 W + (pp 2 IyyW x + W x jT pAa^do^j + p~bu 2 cos 2 AQ,« + ©) 

= pA (vv + p 2 Zp) + Poobu 2 cos 2 AC, — pAg cos <j> (88) 


B. Solution Methods 

Structural dynamic problems can be solved by various computational methods such as finite-element, finite-difference, 
and Galerkin methods. Computational structural dynamic problems are typically formulated in terms of the mass 
matrix, stiffness matrix, and force vector in the form of 19 

[-or (M''+M") - (oG(k)H+ (K s + K c )] U=F (89) 

where M l is an inertial mass matrix, M n is an apparent mass matrix due to noncirculatory forces, K s is a structural 
stiffness matrix, K e is an aeroelastic stiffness matrix due to circulatory forces, F is a force vector, H is a matrix due to 
the phase lag resulting from vortex shedding, and U is a displacement vector. Due to the aeroelasticity, the matrices H 
and K c are non-symmetric. Therefore, the eigenvalues of this system are generally complex. 

Equation (89) is a generalized nonlinear eigenvalue problem due to the presence of G(k) which is a nonlinear 
function of ft). Thus, the eigenvalue solution is an iterative process by first guessing for k, and then solving an 
generalized quadratic eigenvalue problem for ft), which in turn is used to update k. An alternate approach is to simply 
ignore G{k) since G (k) is a small value. The eigenvalue problem then reverts to a linear form which can easily be 
solved. When G(k) is ignored, then the problem is classified as quasi-steady state aeroelasticity. 

There are two types of problems: 1) static aeroelasticity and 2) dynamic aeroelasticity. The static aeroelasticity 
describes physical effects that do not involve dynamic responses of a wing structure such as divergence and control 
reversal. The wing-deflected shape can influence the aerodynamics of a vehicle. The static deflection can be large 
if the wing structure is highly flexible. The static problem can be formulated as a coupled fluid-structure interaction 
problem. The vehicle is modeled by computational fluid dynamics (CFD) method for aerodynamic calculations of 
coefficients and derivatives. The results are used as inputs to the computational structural model to compute the 
vehicle deflected shape. This shape is then used as the new input to the CFD model and the process is repeated until 
the solution converges. 

The dynamic aeroelasticity describes flutter behaviors and transient responses of a wing structure that is subject 
to wind gusts or instantaneous control surface deflections. Dynamic responses of the wing structure can affect the 
overall vehicle control and stability. The solution of the dynamic problem can be solved by implementing the modal 
decomposition method on the eigenvalue solution of the computational structural dynamics. Elastic modes can be 
described by a set of uncoupled, scalar second-order differential equations that can be readily analyzed. 

For a symmetric aircraft configuration for which both wings are identical, two types of elastic modes are present. 
Symmetric modes are those for which the structural deflections of both wings are in the same sense. Anti-symmetric 
modes are those that exhibit structural deflections in an opposite sense between the left and right wings. Figure 8 
illustrates symmetric and anti-symmetric modes. 

Because of the symmetry that exists at the fuselage centerline, only one wing can be analyzed with appropriate 
boundary conditions. 20 The fuselage and tails contribute to the elastic modes as a concentrated mass, half of which is 
located at each of the wing roots. Then for symmetric modes, the boundary conditions at the wing roots must match 
the bending displacement slopes for both left and right wings, and are given by 
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The torsion boundary condition for symmetric modes is null which corresponds to a free-free boundary condition. 
For anti-symmetric modes, the boundary conditions at the wing roots must be zero 
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’ 0 " 

w (0,0 
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_ v ( 0,0 . 
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Anti-Symmetric Modes 



nose-up (down) 


Symmetric Modes 



Fig. 8 - Illustration of Symmetric and Anti-Symmetric Modes 

The symmetric boundary condition is an approximate analysis. An accurate prediction of aeroelasticity requires 
the entire aircraft structure including wings, fuselage, and tails to be modeled. 

The structural deflections obtained from a computational model can be expressed as 
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where the overbar symbol denotes static solutions, ^©^v),/ and ^(©.w.vjj are normalized eigenvectors, qj and rj 
are generalized coordinates for the j — th symmetric and anti-symmetric modes that solve a set of uncoupled scalar 
differential equations 


: 


<ij + c 8 j qj + k 8 j "qj = g S j {u,w,u,w,p,p,8 a ), 
m f' r j + cf" fj + k?* n rj = hf n (u , vv , u , w, p , p , 8 a ) , 


j = 1 
j = 


(93) 

(94) 


The scalar quantities m gen , c gen , and k gen are called generalized mass, damping, and stiffness matrices, respectively, 
and g ge " and h gen are called the generalized forces. They can be computed using the standard procedure in elementary 
vibration theory. 


VIII. Flight Dynamics of Flexible Aircraft 

A flexible aircraft has various elastic modes that can participate in the motion to affect its flight characteristics. 
Wing elastic modes constitute significant structural dynamics of flexible aircraft. In addition, fuselage bending modes 
are also known to affect pitch characteristics. There are other elastic modes such as those due to horizontal stabilizers 
and vertical stabilizer. In totality, all these components contribute to flight characteristics and should be included in 
the equations of motion. The coupled flight-structural dynamics can be quite complicated when all elastic modes are 
accounted for. In this study, we will only focus on the coupling of wing elastic modes and rigid-body flight dynamics. 
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A. Aeroelastic Forces and Moments in Aircraft Reference Frame 


We define the elastic angle of attack as a function of the structural deflections as 
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( 95 ) 


Wing structural deflections affect the lift characteristics of an aircraft. Assuming both the left and right wings are 
of the same geometry, then the elastic angle of attack contributes to the static and dynamic forces and the moments in 
the aircraft reference frame B as 
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where X, Y, and Z are the aircraft axial, side, and normal forces, L, M, and N are the aircraft rolling, pitching, and 
yawing moments, the overbar symbol denotes static forces and moments, and the superscript <J> and l f ; denote dynamic 
forces and moments corresponding to symmetric and anti-symmetric modes, respectively. 

It can be seen that symmetric modes affects forces and moments in the longitudinal direction. On the other hand, 
anti-symmetric modes affects the lateral motion of the aircraft. The aeroelastic forces and moments contribute to 
the flight dynamics of an aircraft in such a way that elastic deflections can adversely affect the rigid-body aircraft 
responses and can result in elastic mode interactions with a flight control system. These interactions necessitate the 
use of aeroelastic mode filters in the flight control design in order to attenuate structural dynamic responses of flexible 
aircraft lifting surfaces. 


B. Propulsive Forces and Moments in Aircraft Reference Frame 

In addition, the elastic modes also affect the propulsive forces and moments generated by wing-mounted engines. The 
propulsive force for the left engine is given by Eq. (72). Assuming a twin-engine aircraft configuration, the propulsive 
force for the right engine is given by 

T R = T [bi - V x b 2 + (WjSinA - ©cosA)b 3 ] ;t=JC£ (99) 

where V, W, and © are defined according to the right wing reference frame C such that V and © are in the opposite 
sense to those deflections in the reference frame D, i.e., V is positive towards the leading edge and © is positive nose 
up. 

The total static and dynamic propulsive forces are then computed to be 

X T = 2 T (100) 

( 101 ) 

X=X E 

where the superscript T S and TA denote thrust forces due to symmetric and anti-symmetric modes, respectively. Thus, 
symmetric modes create a normal force and anti-symmetric modes create a side force due to the combined thrust of 
the two engines. 

The propulsive moment is computed as 

M = (if b, - if b 2 + ifb 3 ) x T L + (if b, + if b 2 + ifb 3 ) x T* (102) 
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where (rf ,rS ,rf ) are the coordinates of the right engine thrust center relative to the aircraft center of gravity in the 
reference frame B. 

Upon evaluation, this yields the total static and dynamic propulsive moments 

M t = 27>f (103) 
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It can be seen that symmetric modes result in an additional pitching moment, whereas anti-symmetric modes create 
both rolling and yawing moments. 

C. Equations of Motion 

The flight dynamic equations with elastic mode and propulsive force interactions can now be written as 
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m{w — qu + pv — gcosdcoscj)) = — CLqS cos a — CoqS sin a cos p + A Z+ (AZj s + AZ 5 ) qj (107) 

7=1 

n 

LxxP T Uv/n' ^x:Pq ~ (Jzz lyy) q r ~Ylyz (*"*" 7~) = CiqSb-\- ^ (A Lj + A L^rj (108) 

7=1 

m 

-IxyP + Iyyq-Iyzr + IyzPq-Ixvqr+(Ixx-Izz)pr + Ixz(p 2 -r 2 ) =C m qSc + M T +XM+Yj (AMJ s + AM*) qj (109) 

7=1 

n 

- IxzP - : V-A + hzT + h z qr - IyzPr + (Jyy ~ hx) pq + Ixy ( q 2 - P 2 ) = C n qSb + £ ANj A r j (110) 

7=1 

where the coefficients with the overbar are for the aircraft, q is the dynamic pressure, c isthe mean aerodynamic chord, 
b is the wing span, S is the reference wing area, and /j is the angle of sideslip. 

These equations constitute six degrees of freedom dynamics that are coupled with structural dynamics equations 
from Eqs. (93) and (94), which are obtained from a computational model. In analyzing the elastic modes, only 
the first few significant modes whose natural frequencies are within a flight control frequency bandwidth are usually 
considered. These elastic modes can become excited by the flight control surface deflections. Examining Eqs. (105) to 
(110) reveals that symmetric modes only affects longitudinal dynamics of aircraft. On the other hand, anti-symmetric 
modes play a role in both lateral and directional dynamics. 

IX. Computation 

In order to establish the equations of motion of flexible aircraft, structural and aeroelastic eigenvalue analysis must 
be performed. Equations of motion of flexible aircraft can conveniently be derived using the Lagrangian equations of 
motion in quasi-coordinates. 1415 This approach allows the derivation of the equations of motion in terms of quasi- 
velocities (time derivatives of quasi-coordinates) in an accelerating reference frame fixed to the aircraft, as opposed 
to the generalized coordinates in a non-accelerating inertial reference frame. The Lagrangian equations of motion 
require the knowledge of three scalar quantities, namely, kinetic energy, potential energy, and the virtual work due to 
the applied forces. In our modeling approach, we regard the aircraft as a flexible multibody system where the aircraft 
components are the fuselage (/), wing (w), horizontal stabilizer (h), and vertical stabilizer (v). To describe the motion 
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of the aircraft, a set of body axes ( x,y,z ) are attached to the undeformed aircraft at a convenient point on the fuselage 
(not necessarily the center of mass of the aircraft), as well as similar axes (, Xj,yi,Zi ), i = w,h,v to the remaining aircraft 
components. The fuselage, right and left wings, right and left horizontal stabilizers, and vertical stabilizer are modeled 
as thin-wall beams, each fixed at its respective root, and subject to flapwise bending displacement W, (x, ) and one 
spanwise generalized torsional displacement ©,■ (x,), i = f. w. h . v. Then the motion of a point on the aircraft can be 
expressed by the rigid-body translation and rotation of the body axes (x ; -,y,-,z,) and by the elastic deformation of the 
flexible bodies relative to their respective body axes. 

The resulting equations of motion include 12 first-order ordinary differential equations for the rigid-body transla- 
tion and rotation, and a partial differential equation for each elastic displacement component. The system is hybrid 
since it includes both ordinary and partial differential equations. The solution method requires a discretization of the 
partial differential equations. To this end, we assume that each elastic displacement can be expressed as a matrix of 
n shape functions multipled by a vector of n generalized coordinates. The shapes functions can be obtained from 
the finite element method, Galerkin method, or any other similar methods. However, for best accuracy at the lowest 
possible value of n, the shape functions used in this paper are the eigenfunctions of a uniform cantilever beam for the 
bending displacements, and the eigenfunctions of a uniform shaft for the torsional displacements. By discretization, 
the partial equations are reduced to a set of first-order ordinary differential equations. 

Aircraft structures generally have a quite complex geometry. For simplicity, the structures of the flexible compo- 
nents are modeled as thin-wall beams with constant thicknesses. To approximate the mass and stiffness distributions 
of the aircraft components, the cross-sectional properties such as geometric center, cross-sectional area, area moments 
of inertia, torsional constant, are computed at some finite number of sections on the respective aircraft components. 
For each flexible body, the geometric center locations are fitted to a straight line which is chosen as the x-axis for 
this body. The z,-axis is chosen to be in the flapwise direction. The aircraft model used in the study is based on the 
NASA Generic Transport Model (GTM). The locations of the geometric centers for the sections and the x,-axes for 
the bodies of this model are shown in Fig. 9. The body axes of the fuselage are treated as the body axes of the aircaft 
and transformation matices are established between the body axes of the aircraft and the body axes of the individual 
aircraft components. Flexural and torsional rigidity distributions of the individual aircraft components with respect to 
their x,-axes are shown in Fig. 10. 



Fig. 9 - Stick Model of Aircraft 

The modal characteristics of the aircraft without aeroelasticity is first estimated by solving the eigenvalue problem 
for the system 


(-© 2 M'' +K s )U = 0 (111) 

where U is the displacement vector for the whole aircraft, including the rigid-body displacements and the elastic 
displacements. 

For aeroelastic eigenvalue analysis, the aircraft is trimmed at a selected flight condition of mach 0.8 and 30,000 ft. 
Then all pertinent aerodynamic derivatives are obtained from an aerodynamic model. This information is used to create 
the aerodynamic stiffness, damping, and mass matrices. To simplify the analysis, the nonlinear effect due to G ( k ) is 
neglected, so that the model is effectively a quasi-steady state model which does not account for the Theodorsen’s 
function due to unsteady aerodynamic lag. 
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Fig. 10 - Structural Rigidity Distributions of Fuselage, Wing, Horizontal Tail, and Vertical Tail 

Table 1 lists the computed structural and aeroelastic modal frequencies and aerodynamic damping associated with 
aeroelastic modes from the eigenvalue analysis. 


Mode 

Structural Frequency, co„ 

Aeroelastic Frequency, co a 

Aerodynamic Damping, 

1 


>> 

II 

o 


2 


A = —0.00165 


3 


A = -1.1423 


4 


0.0088 

0.0138 

5 


0.3681 

0.1508 

6 


0.5498 

0.2868 

7 

7.0633 

6.9637 

0.0120 

8 

7.8161 

7.9244 

0.0019 

9 

7.9600 

8.2228 

0.0044 

10 

10.446 

10.958 

0.0063 

11 

11.029 

11.023 

0.0070 

12 

12.485 

12.445 

0.0075 

13 

13.819 

13.492 

0.0287 

14 

14.285 

13.815 

0.0007 

15 

17.206 

17.218 

0.0134 

16 

33.563 

33.552 

0.0024 

17 

39.272 

39.234 

0.0232 

18 

48.962 

49.277 

0.0099 

19 

49.277 

49.537 

0.0018 

20 

52.126 

53.173 

0.0156 


Since the real parts of the eigenvalues are negative corresponding to positive damping, the aircaft is stable about 
this steady level flight. The first six modes are rigid-body modes, with the first being a purely translation mode of zero 
eigenvalue, the second being the spiral mode, the third being the roll mode, the fourth being the phugoid mode, the 
fifth being the dutch roll mode, and the sixth being the short period mode. The remaining modes are aeroelastic modes 
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with structural deflection involving bending and torsion. The seventh mode is the lowest bending mode and involves 
symmetrical bending of the wing and the horizontal stabilizer about the x — z plane and bending of the fuselage. The 
eighth mode is the lowest torsional mode and involves anti-symmetric bending and torsion of the all elastic members. 
These mode shapes are illustrated in Fig. 3. 



7th mode (bending) at w = 6.96 Hz 8th mode (torsion) at CO = 7.92 Hz 

Fig. 1 1 - First Bending and First Torsion Aeroelastic Modes 

The lowest aeroelastic frequency is about 7 Hz. The aeroelastic stiffness effect does not appear to be too significant 
since the aeroelastic modal frequencies are almost the same as the structural modal frequencies. However, aerodynamic 
damping can have a significant influence in the overalll aircraft stability. 

X. Conclusions 

This paper has presented an integrated flight dynamic modeling method for flexible aircraft. The method combines 
structural dynamics of an equivalent beam model of a flexible wing with rigid-body flight dynamics that accounts for 
coupled effects due to aeroelasticity, inertial forces, and propulsive forces. A formulation of aeroelastic angle of 
attack for the combined chordwise bending, flapwise bending, and torsion is developed that extends Theodorsen’s 
theoretical result. The structural dynamic equations can be solved using the finite-element method to determine static 
and dynamic structural deflections as functions of aircraft states. The elastic modes are decomposed into symmetric 
and anti-symmetric modes with associated generalized coordinates. The standard flight dynamic equations for six 
degree-of-freedom motion then includes the generalized coordinates as additional state variables. These equations 
become coupled with a set of uncoupled second-order differential equations in terms of the generalized coordinates 
that describe the elastic responses of the wing structure. These equations must be solved simultaneously to obtain a 
solution that describes the combined motion of flexible aircraft. A computational aeroelastic model of a generic aircraft 
has been created. Eigenvalue analyses have been conducted to estimate structural modal frequencies and aeroelastic 
frequencies. In the future, using this model, a modal truncation will be performed to generate frequency response 
models that will be coupled with rigid-body flight dynamics. 
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